%% Error Scaling in Numerical Integral Methods


close all

%Limits of integration
A=0;
B=1;

%Function to integrate
%F='sin'
F='exp'

switch F
    case 'exp'
        F=@(x)exp(-x)
        ExactInt=-exp(-B)+exp(-A);
    case 'sin'
        F=@(x)sin(x)
        ExactInt=-cos(-B)+cos(-A);
end

%Limits of integration
A=0;
B=1;

%Range over which to evaulate error
OrdMag=6;

%Setup Arrays
a_b_Array=zeros(OrdMag,1);
Error_MP=zeros(OrdMag,1);
Error_T=zeros(OrdMag,1);
Error_S=zeros(OrdMag,1);

%Find Errors at each sub-interval size
%See modifications to myIntFunc
for nn=1:OrdMag
    N=10^(nn-1); %Number of intervals
    [Nint,a_b]=myIntFunc(F,A,B,'MidPoint',N);
    a_b_Array(nn)=a_b; %sub-interval size (now returned  by 'myIntFunc'
    Error_MP(nn)=abs(Nint-ExactInt);
    [Nint,a_b]=myIntFunc(F,A,B,'Trapezoid',N);
    Error_T(nn)=abs(Nint-ExactInt);
    [Nint,a_b]=myIntFunc(F,A,B,'Simpson',N);
    Error_S(nn)=abs(Nint-ExactInt);
end


%Look at error vs sub-interval size
figure
loglog(a_b_Array,Error_MP)
hold on
loglog(a_b_Array,Error_T)
loglog(a_b_Array,Error_S)
loglog(a_b_Array,a_b_Array.^2,'--') %(a-b)^2 reference
loglog(a_b_Array,a_b_Array.^4,'--') %(a-b)^4 reference
loglog(a_b_Array,ones(size(a_b_Array))*eps(ExactInt),'--')
xlabel('Interval Size |a-b|')
ylabel('Error from Exact')
legend('MidPoint','Trapezoid','Simpson','(a-b)^2','(a-b)^4','eps')
title(char(F))
axis equal

%Note that we start to get funny behavior when the error is reaches 'eps'
%This is coming from the resolution of floating point numbers
%doc eps

%% Some examples of machine precision limits

format long

%The next largest number than 1 which can be represented by 'double'
%floating point numbers is
1+eps(1)

(1+eps(1))-1

%We have problems when we add numbers, or numbers with significant digits,
%below eps
(1+eps(1)/2)-1

1.111111111*10^-14
X= 1+1.111111111*10^-14
1-X

%% A demonstration of the Central Limit Theorem

%web('https://en.wikipedia.org/wiki/Central_limit_theorem#Classical_CLT')

%Change NSums from 1 to a large number
NSums=10
NRepeats=100000;

%Random numbers from a uniform distribution on -1 to 1
Samples=(rand(NRepeats,NSums)-.5)*2;

%Add Random Numbers together
S=sum(Samples,2);

%Note as NSums -> Inf, the distribution is more Gaussian
figure
hist(S,100)
xlabel('S')
ylabel('Counts')

%% A demonstration of the Central Limit Theorem with a more extreme distribution

%Change NSums from 1 to a large number
NSums=1000
NRepeats=10000;

%Each number is either -1 or 1 
Samples=round(rand(NRepeats,NSums))*2-1;

S=sum(Samples,2);

figure
hist(S,(min(S)-1:max(S)+1))
xlabel('S')
ylabel('Counts')


%% Random Walk

% What we have done above is the same as simulating the end positions of 
% random walks

% web('https://en.wikipedia.org/wiki/Random_walk')

%Let's look at a 2D random walk vs time step

NSums=1000000; %Lots of steps
NRepeats=1; % only one walker

SamplesX=round(rand(NRepeats,NSums))*2-1;
SamplesY=round(rand(NRepeats,NSums))*2-1;

%Calculate the displacement vs time step
%help cumsum
CumSX=cumsum(SamplesX,2);
CumSY=cumsum(SamplesY,2);

figure
plot(CumSX(1,:),CumSY(1,:))
xlabel('X')
ylabel('Y')

%Try running this a few times and look at the differing results





